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The phase diffusion in a self-sustained oscillator, which produces oscillator's spectral linewidth, 
is inherently governed by a nonlinear Langevin equation. Over past 40 years, the equation has been 
treated with linear approximation, rendering the nonlinearity's effects unknown. Here we solve the 
nonlinear Langevin equation using the perturbation method borrowed from quantum mechanics, and 
reveal the physics of the nonlinearity: slower phase diffusion (linewidth narrowing) and a surprising 
oscillation frequency shift that formally corresponds to the Lamb shift in quantum electrodynamics. 

PACS numbers: 02.50.Ey, 31.30.J- 
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PROBLEM STATEMENT 

Self-sustained oscillators are a major research subject 
in the broad context of science and technology. Elec- 
tronic oscillators [H Q, masers/lasers SB, neural cir- 
cuits [H*] , and human circadian rhythms @ are examples 
from electronics, optics, and physiology. Such oscillators 
sustain oscillations on limit cycles by compensating par- 
asitic energy loss [l!]. Noise disturbs oscillation, causing 
amplitude and phase errors. The amplitude error that 
puts oscillation off the limit cycle is constantly corrected 
by oscillator's tendency to return to its limit cycle. By 
contrast, the phase error on the limit cycle accumulates 
without bound, as no mechanism to reset phase exists. 
Thus phase undergoes Brownian motion, diffusing along 
the limit cycle. The oscillator signal with frequency wq 
may then be written as v{t) — vq cos[u!ot + (j){t)] ignoring 
the amplitude error, (jj^t) is the phase diffusion. Due to 
(j){t), w(t)'s spectrum is broadened about wg. The phase 
diffusion and spectral broadening are among the most 
essential aspects of oscillators' dynamics and quality. 

The phase diffusion due to the disturbance by a Gaus- 
sian white noise G{t) is governed by an inherently non- 
linear Langevin equation 0, @| '■ 



d(l)/dt = cos(wot + <l>)G{t). 



(1) 



This is the standard Langevin equation for Brownian 
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motion, except the periodic modulation cos(u;oi -I- <j>) 
that depends on the oscillator's state, ujot + (j>{t). This 
state-dependent modulation is a hallmark property of the 
phase diffusion 0, S]- it arises, as the phase error for a 
given disturbance varies with where the disturbance oc- 
curs on the limit cycle. This is shown in Fig. [T]for an LC 
oscillator. During the time interval At, current i in the 
LC tank is disturbed by noise GAt produced by parasitic 
resistor R. The resulting phase error A<j) is a projection 
of GAt along the tangential direction of the limit cycle, 
which depends on the limit cycle position. The modula- 
tion can be a general wg-periodic function, for the noise 
intensity can also vary with the oscillator's state and the 
limit cycle can assume a non-circular shape, but the co- 
sine modulation captures the essence. 

When (p diffuses according to ([T]), probability den- 
sity p{4>, t) evolves according to the Fokker-Planck equa- 
tion, dp/dt = {d^ /d(f>^)[{a^/ 2) cos^ {ujot + (f>)p]- Here 
is G{tys power spectral density: {G{ti)G{t2)) = 
(T"(5(ii — ^2), where (•) is ensemble average. The instanta- 
neous phase diffusion rate, cr^ cos^(a;oi -I- 0)/2, identified 
in the Fokker-Planck equation, depends on the oscilla- 
tor state. Its time average, after ignoring the nonlinear 
^-dependence, is I? = (7^/4, which gives a sense of the 
average diffusion rate (the exact value can differ from D 
due to the 0-dependence, as seen later). D and wq are 
oscillator's two characteristic frequencies. 

Due to the nonlinear 0-dependence, ([1]) has no closed- 
form solution. In 1960s, phase diffusion was studied in 
the low-noise regime, D <^ wn, to which electronic os- 
cillators and lasers belong (iHil. These works dropped 
the nonlinear ^-dependency in ([T]), which allowed time- 
averaging of the remaining modulation to reduce ^ to 



d(f>/dtKil/V2)G{t). 



(2) 



FIG. 1: (a) Circuit diagram of a self-sustained LC oscillator, 
(b) Limit cycle and state-dependent phase error. 



This linear approximation works, as for D <^ ujq many 
cycles of oscillation occur before phase diffuses apprecia- 
bly. The corresponding simplified Fokker-Planck equa- 
tion gives a normally distributed p{4>^t) with (0^(t)) — 
2Dt. This diffusion with rate D leads to the well-known 
Lorentzian spectrum with half-width D at half-power Q . 
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For over 40 years, however, the general nonhnear 
Langevin equation ([1]) has escaped solution, and a cou- 
ple of fundamental questions on phase diffusion have re- 
mained unanswered. First, while the simplified model 
([2]) is valid in the low-noise regime D ujq (electronic 
oscillators, lasers) 0, the lack of solution to the gen- 
eral equation ([1]) has obscured when/how the simplified 
model fails, shadowing confidence of electronics and laser 
community in the simplified model. Second, in the high- 
noise regime D ^ cjo where the nonlinear 0-dependency 
cannot be ignored, what are the physical effects of the 
nonlinearity? Not only is this in itself a fundamental 
question on phase diffusion, but the answer can be useful 
in studying timing accuracy of low-frequency high-noise 
oscillators, e.g., neural circuits circadian rhythms Q. 

Here we address this long-standing problem by solv- 
ing the nonlinear Langevin dynamics (jlj (Sec. 2). We 
thus answer both questions above (Sec. 3), drawing the 
boundary of the linear model and revealing the physics of 
the nonlinearity: slower phase diffusion (narrower spec- 
trum) and a surprising oscillation frequency shift. Also 
surprising is the formal correspondence of the frequency 
shift in classical oscillators to the Lamb shift in quantum 
electrodynamics (Sec. 3); this finding highlights, with a 
widely-used physical system as a concrete example, the 
intimate link between quantum and stochastic systems 
known from works by Nelson and Feynman 11 1 . 



SOLUTION VIA FEYNMAN-KAC EQUATION 

The power spectral density of v(t) = vq cos{ujQt + (/)),& 
measurable quantity capturing phase diffusion, is Fourier 
transform of the autocorrelation Rv{t) — (^v{t)v{t + r)). 
Thus our task is to compute . Our approach is to relate 
Ry to a function (which we call /) governed by Feynman- 
Kac equation, and solve it using the perturbation method 
borrowed from quantum mechanics. 

We deal with the total phase ijj{t) = uj^t -t- (i){t). The 
nonlinear Langevin equation for '0(i)) matched to ([TJ, is: 

d'0M = wo + cos-0-G'(i); (3) 

and the autocorrelation with v{t) — vq cos ip{t) reads 

Ry{T) = it;2Re[(e'W'(*+-)-'^(*)l) + (e^['>(*+-)+'^(*)l)]. (4) 

Here (e*['/'(t+^)T0(t)]) ^ (^^T^m(^e''f'^t+r)^^^^■^^^^^^^ 

conditional probability. Let /(a:,r) = (e*'^(*+^) |i/'(i) = 
x) be the ensemble average of e"'''^*+'^\ given ■0(t) = x e 
[0, 27r], and p{x) be the stationary probability density of 
Tpit) — X for t large enough. Then we have 

l^^mt+r)^W)\ ) = T)eT"p(x)dx. (5) 

We have reduced computation of R^ to that of /(x, t) 
and p(x). Compared to /, p is easier to calculate (p has 
a closed-form expression |12|). and has no major physical 



effects as seen later, so we focus on /. / decays with r, 
for an ensemble of oscillators with the same initial phase 
X dephase over time due to phase diffusion (decay rate = 
dephase/diffusion rate), thus, / captures phase diffusion 
and is physically meaningful. Mathematically, given ([3]), 
/ satisfies the following Feynman-Kac equation for t > 0: 



dr 
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dx 
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D cos 2x 



dx'^ 
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(6) 
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This does not contain t, which is why t-dependency was 
dropped in / and R-^. The remaining job is to solve 
It has no closed-form solution due to the H' term, whose 
origin is the nonlinear i/i-dependency in the Langevin 
equation ([3]). As ([6]) has the same form as Schrodinger 
equation, we solve ([6]) using the perturbation method of 
quantum mechanics, separating it into unperturbed (Hq) 
and perturbed {H') terms. When D <C ujq, H' is negligi- 
ble, corresponding to the simplified model 

Let e„ and g„ (x) be the n-th eigenvalue and eigenstate 
of®: {Hq + H')gn = enQn- Quantization n = ±1, ±2, ... 
is due to periodic boundary condition (7„(0)=(7„(27r). e„ 
and gn can be obtained via the perturbation technique. 
As gn evolves as e'^"^5„, /(x,r) = c„e'^"^g„(x). Co- 
efficients c„ are set by initial condition /(x, 0) ~ e". 

The unperturbed equation H^gn — £ngn with the peri- 
odic boundary condition yields the Oth-order g„ and e„: 



.(0) 



inujQ - 



(7) 



Each eigenstate is a harmonic mode of oscillation; its fre- 
quency nwo and decay rate n^D are lm[el'^^] and Re[ei°'']. 
In this unperturbed case, initial condition /(x, 0) = is 
an eigenstate, gf^\x), of Hq, so f{x,T) = e*^! ^'^g^i\x) ~ 
q-Dt^i(ujot+x) -^^rith no other eigenstates excited. 

For higher-order solutions, we obtain matrix elements 
of H' with the unperturbed eigenstates gn \x) as basis: 



H' 



i^Cgii:^*H'gl:^Ux^--in^D/2)S^,n 



±2- 



(8) 



This off-diagonal matrix alters eigenstates g„ and eigen- 
values e„ of Hq + H' from the unperturbed ones, whose 
calculation to the lowest order we describe now. In the 
lowest-order correction, gn^ is coupled only with g^^2 



given dH), thus, g„ = 5!°' + a„,„+25^+2 + "■n,n-2g, 



(0) 



.(0) 



with an,n±2 — -ff'i±2.n/(^«'' ~ '^n±2) obtained from the 
perturbation formula of quantum mechanics. Also e„, 
or the frequency and decay rate, is altered due to gi^^'s 



(0) 



coupling with g^^j.2 in its lowest-order correction: 



H' H' 

E mn nm 
JO) _ (0) 
m=n±2 "^n tm 



With gn and e„ available to the lowest order, we can 
find /(x, r) to the same order by using the initial con- 
dition, /(x,0) = e'^. As = ^^^■'(x) is not an eigen- 
state of Hq + H' , it can be expanded in terms of g^s, 
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the eigenstates of Hq + H'; to the lowest order, we find 
e'^ = cigi + c_ig_i + C353 with ci = 1, c_i = D/ (iAujo), 
C3 = D/{—i4ujQ + 16D) (|ci| is far larger than |c_i| and 

jcal for D <C wq, as f(x, 0) — g^^\x) is still close to 
gi{x) for D <^ ljq) [l^l- We can then readily write 
f{x,T) = cie"i^gi(x) + c_ie'-i'^si_i(a;) + Cse'^^^gsix) to 
the lowest order with no other eigenstates excited. The 
e'^±i'^ terms dominate as they have slowest decay rate 
(~ ri^D = D). Ignoring the e'^^'^ term with far higher 
decay rate ri^D = 9D) and using e„ = el„, ffn = 



fix, r) « cie'i^5i(a;) + c^ie"'^^ gl{x). 



(10) 



We can extend the calculation to any higher order [12] ■ 
In higher orders, / contains more excited eigenstates, but 
(1101) still holds (with more accurate ei, gi, and c±i), given 
the slowest decay of the two terms (secondarily, |ci| is 
larger than other |c„|'s due to the initial condition). 

Using ([To]) in ([5]) yields Rv To the lowest order, it is: 



Rvir) 



1^ (r,-a)|r| 



Here a = Re[ei — ef''] and /? = Im[ei — e^'^''] become the 
lowest order correction of ei according to ^ : 

a = {9/2)[D^/{mD^ + ujI)]D; (12a) 
f3 = {D/8)[9Dujo/{l6D'^ + c^l) ~ D/ua], (12b) 

which represent decay rate reduction and frequency shift, 
caused by g^'''s couphng with g'^^l, g^\ For ([TT1) . we 
approximate p{x) of ([5]) as uniform, which is accurate for 
D <^ loq. The exact p{x), found in closed form by solving 
the stationary Fokker-Planck equation corresponding to 
([3]), slightly alters coefficients of the terms in while 
the rest, including the physically important a and /3 that 
are from /(x, r), remain the same 1^. In higher orders, 
the functional form of ([Tl]) is maintained, as (jlOp is still 
used, while the coefficients, a, and /3 are more accurate 
due to higher-order corrections of ei and gi. 

SLOW DIFFUSION AND FREQUENCY SHIFT 

The autocorrelation Ry in ([Tl]) captures physical con- 
sequences of the nonlinear Langevin dynamics ([T]). First, 
the decay rate D — a, or the linewidth of the power spec- 
trum corresponding to R^, is the phase diffusion (de- 
phase) rate, which is smaller than the rate D obtained 
from the approximate linear model Second, the os- 
cillation frequency loq + /3 seen from (jlip is surprisingly 
shifted from the natural oscillation frequency wq by /3. In 
the low- noise regime [D ^ luq) where the lowest-order a 
and /3 in (fT2|) are highly accurate, a/D a (D/ujo)'^ <^ 1 
and P/ujQ (X (Z?/wo)^ ^ 1, thus the true nonlinear dy- 
namics dD) and simplified model ^ quadratically con- 
verge as D decreases. In contrast, if D increases towards 
Wo, a and (3 grows appreciable: for D = ujq, a/D = 0.26 



cos(a;o + P)t ~ ^ sin(a;o + /3)|t| 

(11) 



from (|12ap . and a/D = 0.31 by the n-th order calculation 
until convergence; the nth-order /3 is 8% of wq. The anal- 
ysis decisively justifies the linear model ^ for _D ^ wq, 
finding its boundary set by the quadratic convergence, 
thus, enhancing confidence in the linear model used in 
lasers and electronic oscillators. It has been also found 
that the physical effects of the nonlinearity are the slower 
diffusion and frequency shift, and they are conspicuous 
in the high-noise regime {D ^ ujq). 

The slower phase diffusion may be understood as fol- 
lows. The instantaneous diffusion rate 2Dcos^{ujQt + (p) 
is zero when ~ — wpi + n7r/2 with odd integer n. These 
0-points, which change with time, act as barriers that 
prevent phase diffusion across them. When is small, 
the barrier points tend to move slowly or be more fix- 
ated, making their effects more potent. In the extreme 
case of Wo = 0, cannot diffuse to the region {tt/2, 37r/2) 
starting from (t){t = 0) G [— 7r/2, 7r/2]. Overall, this bar- 
rier action leads to slower phase diffusion. 

The frequency shift, no matter how small, is surpris- 
ing and interesting from the perspective of time-domain 
dynamics, as noise disturbs phase in positive and nega- 
tive directions with equal probability and no such shift 
is expected. In our calculation using eigenstates, the fre- 
quency shift comes out as part of the eigenvalue shift. It 
arises from the interactions among the harmonic modes 
{gl?) caused by the Z? cos 2a; term {H') in ([6]), or the 
matrix element as seen in ([9]). The calculation re- 

veals that both the noise fluctuations captured by D and 
nonlinear modulation captured by cos 2x are the origin 
of i?,^„, hence, indispensable for the frequency shift. 

Also interesting is the formal correspondence of this 
frequency shift in classical oscillators to the Lamb shift 
in quantum electrodynamics. In our frequency shift, 
noise fluctuations {D) and the nonlinear modulation 
(cos2j:) play the role of vacuum fluctuations (A) and 
the non-diagonal dipole operator (p) in the Lamb shift. 
The latter two enter the interaction Hamiltonian H' = 
— (e/2mc)(p • A + A • p), which couples different atomic 
states to produce the Lamb shift jl3| : 



E 



m^n photon 



E„ 



hjjj 



(13) 



The summation is over both atomic states and virtual 
photon (vacuum fluctuation) states. Notice its similarity 
to ([9]), except the missing haj term in ([9]). This is sensi- 



(b) y (£• = 0) 





FIG. 2: Feynman diagrams of the lowest-order process for (a) 
the Lamb shift and (b) our frequency shift. 
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FIG. 3: (a),(b) Numerical (e'[*('+"'-'''(*)l) for D = a;o; 
(c),(d) Numerical a, /3 (cross marks) vs. analytical a, (3 
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FIG. 4: Power spectra from Monte-Carlo simulations of U]) 
(solid) and @ (dashed), ojq = 1: (a) D = ojo; (b) D = a;o/4. 

ble, as the "virtual photons" in our case has zero energy, 
for phase disturbance does not incur energy exchange. 
Figure [5] illustrates the correspondence using Feynman 
diagrams of the lovifest-order process. 

Atomic level shifts per se are common in quantum me- 
chanics, but vi^hat makes Lamb shift unique is its origi- 
nation from fluctuations (our focus in Lamb shift is its 
physical origin rather than its computation using renor- 
malization and relativistic corrections). This is why our 
frequency shift, which also occurs due to fluctuations, 
uniquely corresponds to Lamb shift. This finding can be 
put in the fundamental context of the intimate formal 
link between stochastic and quantum dynamics, estab- 
lished by Nelson [l^], Feynman the link arises as 
both deal with the evolution of probability-related func- 
tions. In this context, it is interesting to seek for a real 
physical system as a stochastic analogue of Lamb shift. 
The phase diffusion with the unique natural nonlinearity 
in classical oscillators offers one such stochastic analogue. 



NUMERICAL CALCULATIONS 

To verify our solution to the Feynman-Kac equation, 
we solve ^ using the finite difference method, and cal- 
culate (^e'^'4'it+r)~i'{t)]'^ of ^ The nearly linear relations 

of log|(e*['^(*+^)-'^(*)l)| and arg(e*['''(*+^)-'^(*)l) to r for 
Z? = Wq in Figs.|3Ua),(b) confirm that a single exponential 
g~(D-a)Tgi(wo-i-/3)T jj^figgfi dominates the autocorrelation. 
The location of log |(e*['^(*+^)-'^(*)l)| above the line of 
— Z3t, and the positive slope of arg<^e't^'-*+'^^~'^'^*-*l) —ujqt, 
confirm the reduced diffusion rate and frequency shift. 
Figures [3]Jc),(d) plot a and /3 from our analytical calcu- 
lations to the lowest order ([T^ and the nth order (up to 
convergence at rt ~ 10 [Hj), against their numerically ob- 
tained values. Our analysis agrees well with the numeri- 
cal results; the discrepancy between (|12p and the numer- 
ical result at large D is readily corrected by higher-order 
calculations. Note also that a/D and fi/cuo diminish as 
{D/uj[)Y with decreasing D as shown in Figs. [3][c),(d). 

For further confirmation, we Monte-Carlo simulate the 
nonlinear Langevin dynamics ([T]). The resulting power 
spectra Sv{i^) (Fig. [4|) are sharper/narrower than the 
simplified model, confirming the slower phase diffusion: 
the rate reduction is 31% in otir analysis and 27% in the 
Monte-Carlo simulation for D — uif)] 11% and 12% for 
D = wo/4. The spectrum peak exhibits a definitive blue 
shift in centre frequency (Fig. HJ . The shift for D = ojq/A 
is 2% in Monte-Carlo simulation, and 3% in our analysis. 
For D = ujQ, the spectrum peak does not exactly corre- 
spond to ujQ + /3 due to the large D, thus, reading the 
spectrum requires caution, but the extracted /?, which 
is a good representation of the frequency shift, matches 
well between the simulation (8%) and our analysis (8%). 
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